******************************************************************************
*										
*                              Boese & Eberhardt							
*            
*                             Facets of Democracy
*													
******************************************************************************
*                   PCDID - heterogeneous model (1959-2018)
*			    (Specification with exports/trade & pop growth)
* 				   High and mid-level indices: (1) Lib_dem,
*			       (2) Polyarchy, (3) Liberal component
*                   Varying the cutoff: Alpha and Chi2 tests
******************************************************************************
*                       Using standardized thresholds 
******************************************************************************
*                           Created: 5th June 2024					
*                      Last Changed: 5th June 2024					
******************************************************************************
*                              Markus Eberhardt						
******************************************************************************
*                 School of Economics, University of Nottingham,			
*                  University Park, Nottingham NG7 2RD, England			
*                     email: markus.eberhardt@nottingham.ac.uk				
*                   web: https://sites.google.com/site/medevecon/			
******************************************************************************

clear all
set matsize 11000

***
*** Paths
***

global path "D:/Dropbox/Facets of Democracy"
global path "/Users/lezme/Dropbox/Facets of Democracy"
global rawpath "/Users/lezme/Dropbox"

global data "$path/Data"
global dofolder "$path/Do_files"
global otherdata "$path/Data"
global outputfolder "$path/Output"
global texfolder "$rawpath/Apps/Overleaf/Facets of Democracy"

global xtmgpath = "D:/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"
global xtmgpath = "/Users/lezme/Dropbox/Literature/econometric issues/Nonstationary Panel Metrics/Stata"

****
**** Get merged dataset and prepare variables
****

use "$data/madd_ert_vdem_2021.dta", clear
tsset nwbcode year

* Maddison: per capita GDP (in logs)*100 and population growth rate
gen y = ln(gdppc)*100
gen lpop = ln(pop)
gen dlpop = 100*d.lpop

* DOTS data: export/trade
gen ex = ex_trade*100

* Exclude years prior to 1959
xtreg y dlpop ex v2x_libdem if year>=1949, fe
* n=7,643, N=157, avg T=48.7 (min 12, max 60); 1959-2018

drop csum 
gen c=1 if e(sample)
by nwbcode: egen csum=sum(c) if c==1
gen sample= (e(sample))

sum v2x_libdem  v2x_polyarchy v2x_liberal if sample==1

* V-Dem: democracy dummies from high and mid-level indicators
local z "libdem polyarchy liberal"
foreach k of local z{
	sum v2x_`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st_`k' = (v2x_`k'-mean`k')/sd`k' if sample==1
}

* Positive values
local z "0 25 50"
foreach k of local z{
	local norm_`k' = `k'/100
	gen dem_ld_`k' = (st_libdem>=`norm_`k'') & !missing(v2x_libdem)
	label var dem_ld_`k' "Libdem Threshold > mean + .`k' sd"
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	label var dem_poly_`k' "Polyarchy Threshold > mean + .`k' sd"
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal) 
	label var dem_lib_`k' "Liberal Component Threshold > mean + .`k' sd"
}

* Negative values
local z "50 25"
foreach k of local z{
	local norm_neg`k' = -`k'/100
	gen dem_ld_neg`k' = (st_libdem>=`norm_neg`k'') & !missing(st_libdem)
	label var dem_ld_neg`k' "Libdem Threshold > mean - .`k' sd"
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_neg`k'') & !missing(st_polyarchy)
	label var dem_poly_neg`k' "Polyarchy Threshold > mean - .`k' sd"
	gen dem_lib_neg`k' = (st_liberal>=`norm_neg`k'') & !missing(st_liberal) 
	label var dem_lib_neg`k' "Liberal Component Threshold > mean - .`k' sd"
}


* Positive values
local z "125"
foreach k of local z{
	local norm_`k' = `k'/1000
	gen dem_ld_`k' = (st_libdem>=`norm_`k'') & !missing(v2x_libdem)
	label var dem_ld_`k' "Libdem Threshold > mean + .`k' sd"
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	label var dem_poly_`k' "Polyarchy Threshold > mean + .`k' sd"
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal) 
	label var dem_lib_`k' "Liberal Component Threshold > mean + .`k' sd"
}


* Negative values
local z "125"
foreach k of local z{
	local norm_neg`k' = -`k'/1000
	gen dem_ld_neg`k' = (st_libdem>=`norm_neg`k'') & !missing(st_libdem)
	label var dem_ld_neg`k' "Libdem Threshold > mean - .`k' sd"
	gen dem_poly_neg`k' = (st_polyarchy>=`norm_neg`k'') & !missing(st_polyarchy)
	label var dem_poly_neg`k' "Polyarchy Threshold > mean - .`k' sd"
	gen dem_lib_neg`k' = (st_liberal>=`norm_neg`k'') & !missing(st_liberal) 
	label var dem_lib_neg`k' "Liberal Component Threshold > mean - .`k' sd"
}


*** RANGES OF THE ROBUSTNESS EXERCISE

di "LIBDEM: -1/4: "	 -.25*.2819623 "  -1/8: " .3544331-.125*.2819623 "  mean: " .3544331 "  +1/8: " .3544331+.125*.2819623 "   +1/4: ".3544331+.25*.2819623
**
*LIBDEM:	-1/4:	.28394253	-1/8:	.31918781	mean:	.3544331	+1/8:	.38967839	+1/4:	.42492368
**

di "POLY: -1/4: ".4514585-.25*.2889042 "  -1/8: " .4514585-.125*.2889042 "  mean: " .4514585 "  +1/8: " .4514585+.125*.2889042 "   +1/4: ".4514585+.25*.2889042
**
*POLY: -1/4: .37923245  -1/8: .41534547  mean: .4514585  +1/8: .48757153   +1/4: .52368455
**

di "LIB: -1/4: " .5545428-.25* .2894293 "  -1/8: " .5545428-.125*.2894293 "  mean: " .5545428 "  +1/8: " .5545428+.125*.2894293 "   +1/4: ".5545428+.25*.2894293
**
*LIB: -1/4: .48218547  -1/8: .51836414  mean: .5545428  +1/8: .59072146   +1/4: .62690013
**


* ROW: democracy dummy
gen dem_row2 = (v2x_regime>=2) & !missing(v2x_regime)
label var dem_row2 "ROW - Regime Change 1 to 2"

* PolityIV: democracy dummy
gen dem_p4_1 = (polity2>=1) & !missing(polity2)
gen dem_p4_6 = (polity2>=6) & !missing(polity2)

* Boix
label var dem_boix "Boix et al Democracy indicator (1959-2014)"

****
**** PCDID 
**** 

estimates clear 

local i = 1

**** Select the variables, sample period, factor-augmentation 
global xvars1 		"ex dlpop"	// X variables used to construct factors 
global xvars2 		"ex dlpop"	// Xs included in the PCDID (leave blank for 'plain vanilla' model)
global factors 		"6"			// max # of factors 
global uhat 		"uhatMG"	// uhatMG: residuals from heterog regression, uhatP: pooled FE
global demo_list	"p4_1 p4_6 row2 ld_neg25 ld_neg125 ld_0 ld_125 ld_25 boix poly_neg25 poly_neg125 poly_0 poly_125 poly_25 lib_neg25 lib_neg125 lib_0 lib_125 lib_25" 
global demo_list2	"dem_p4_1 dem_p4_6 dem_row2 dem_ld_neg25 dem_ld_neg125 dem_ld_0 dem_ld_125 dem_ld_25 dem_boix dem_poly_neg25 dem_poly_neg125 dem_poly_0 dem_poly_125 dem_poly_25 dem_lib_neg25 dem_lib_neg125 dem_lib_0  dem_lib_125 dem_lib_25" 
 
run "$xtmgpath/xtmg.ado"

* Start of loop for different end years 
local zzz "2018"
foreach kkk of local zzz{
	global endyear 		"`kkk'"			// 2018

	* Start of loop for different start years 
	local zz "1959"
	foreach kk of local zz{
		global startyear 	"`kk'"		// 1949 

		* Start of loop for different democracy indicators 
		local z "$demo_list"
		foreach k of local z{
			global democracy 	"`k'"		// pick democracy indicator (V-Dem indices)

			* Start of factor loop: how many factors to be included in the PCDID
			forvalues f=1/$factors{
				preserve
				quietly{
					tsset nwbcode year
					keep if year>=$startyear & year<=$endyear
					gen dem_dum = dem_$democracy
					sort nwbcode year
					* Cross the threshold
					gen count_$democracy = 1 if ((dem_dum==0 & l.dem_dum==1) | (dem_dum==1 & l.dem_dum==0))
					by nwbcode: egen dem_dum_count_$democracy = sum(count_$democracy)
					replace count_$democracy = 0 if (dem_dum==0 & l.dem_dum==1)
					* Cross the threshold from below
					by nwbcode: egen trudem_dum_count_$democracy = sum(count_$democracy)
					by nwbcode: egen dem_dum_exposure_$democracy = sum(dem_dum)
					by nwbcode: egen sample_obs_$democracy = sum(sample)
					gen autoc_exposure_$democracy = sample_obs_$democracy - dem_dum_exposure_$democracy
					xtlogit dem_dum y if sample==1, fe iterate(1)
					* Treat are those which experienced crossing the threshold
					gen treat = (e(sample))
					replace treat = . if sample==0
					* Below are those which are always below the threshold
					gen below = (!e(sample) & sample==1 & !missing(dem_$democracy))
					replace below = . if sample==0
					* Always above threshold (adjust)
					gen above = (below==1 & dem_$democracy ==1)
					replace above = . if sample==0
					* Always below threshold (adjust)
					replace below = . if below==1 & dem_$democracy==1
					
					*** Estimate linear heterog regression of y on x in the control group and collect residuals
					xtmg y $xvars1 if below==1 & sample==1 & above==0, res(uhatMG)
					local below = e(N_g)
					local belowobs = e(N)
					*** Estimate linear pooled regression of y on x in the control group and collect residuals
					xtreg y $xvars1 if below==1 & sample==1, fe
					predict uhatP, e
					*** Compute residual period average in control group 
					sort year 
					by year: egen uCT = mean($uhat)
					sort nwbcode year					
					*** Extract factors from the residuals - decide which version ($uhat) to use
					xtreg $uhat if sample==1, fe
					gen t_max = e(g_max)
					gen N_max = e(N_g)
					if t_max > N_max{
						regife $uhat if sample==1, f(cfactor = year floading = nwbcode, `f') maxiterations(1)
					}
					else if  N_max > t_max {
						regife $uhat if sample==1, f(floading = nwbcode cfactor = year, `f') maxiterations(1)
					}
					drop t_max N_max
					
					*** Create factors for all countries
					forvalues l=1/`f'{
						sort year
						by year: egen cfactor`l'N = mean(cfactor`l')		
						drop cfactor`l' 
					}
					drop floading*
					
					*** Test of controls ***
					xtlogit dem_$democracy y $xvars2 ///
						if sample==1 & treat==1, fe iterate(1)
					xtmg dem_$democracy $xvars2 cfactor*N ///
						if e(sample)
					testparm $xvars2
					local chi2_stat = r(chi2)
					local chi2_p = r(p)

					*** Alpha test auxiliary regressions ***
					xtmg y uCT ///
						dem_$democracy ///
						$xvars2 ///
						if sample==1 & treat==1
					mat alphas = e(betas)
					svmat alphas, name(alphas)
					gen alpha01 = alphas1-1
					replace alpha01 = . if (alphas2)==0
					mkmat alpha01, mat(alpha01)
					set seed 280115
					reg alpha01
					mat a1 = e(b)
					mat vv1 = e(V)
					scalar alphaest1 = a1[1,1]
					local alpha_est1 = alphaest1
					scalar vv1s = vv1[1,1]
					scalar se1s = sqrt(vv1s)
					scalar alphat1 = alphaest1/se1s
					local alpha_tstat1 = alphat1
					
					*** Run PCDID model ***
					eststo: xtmg y dem_$democracy $xvars2 ///
						cfactor*N ///
						if sample==1 & treat==1, res(resid) robust
					mat betas = e(betas)
					estadd scalar treated = e(N_g)
					estadd scalar obs = e(N)
					local treated = e(N_g)
					local obs = e(N)
					mat b1 = e(b)
					mat v1 = e(V)
					scalar b_1 = b1[1,1]
					scalar se_1 = sqrt(v1[1,1])	
					scalar mse = e(sigma)
					estadd scalar mse = e(sigma)
					capture xtcd resid if e(sample)
					capture scalar cdtest = r(pesaran)
					capture estadd scalar cdtest = r(pesaran)
					estadd scalar fac = `f'
					estadd scalar x_test_chi2 = `chi2_stat'
					estadd scalar x_test_p = `chi2_p'
					estadd scalar alpha_b1 = `alpha_est1'
					estadd scalar alpha_t1 = `alpha_tstat1'
					estadd scalar controlN = `below'
					estadd scalar controlobs = `belowobs'
					estadd scalar final_year = `kkk'
					estadd scalar start_year = `kk'
					est store model`i'
					
					keep if e(sample) & treat==1
					tabstat nwbcode if e(sample) & treat==1, by(wbcode) save
					tabstatmat nwbcodestat, nototal
					tabstat year if e(sample) & treat==1, by(wbcode) stats(min max) save
					tabstatmat yearstat, nototal
					tabstat dem_dum_count_$democracy trudem_dum_count_$democracy dem_dum_exposure_$democracy ///
						sample_obs_$democracy autoc_exposure_$democracy, by(wbcode) stats(mean) save
					tabstatmat exposure, nototal
					drop _all
					svmat nwbcodestat, name(nwbcode)
					rename nwbcode1 nwbcode
					svmat yearstat, name(year)
					rename year1 startyear
					rename year2 endyear
					svmat exposure, name(count)
					rename count1 flips
					rename count2 democratisations
					rename count3 exposure
					rename count4 total_obs
					rename count5 nonexposure
					svmat betas, name(beta)
					gen wbcode=""
					replace wbcode="AFG" if nwbcode==	2
					replace wbcode="AGO" if nwbcode==	3
					replace wbcode="ALB" if nwbcode==	5
					replace wbcode="ARE" if nwbcode==	7
					replace wbcode="ARG" if nwbcode==	8
					replace wbcode="ARM" if nwbcode==	9
					replace wbcode="AUS" if nwbcode==	12
					replace wbcode="AUT" if nwbcode==	13
					replace wbcode="AZE" if nwbcode==	14
					replace wbcode="BDI" if nwbcode==	15
					replace wbcode="BEL" if nwbcode==	17
					replace wbcode="BEN" if nwbcode==	18
					replace wbcode="BFA" if nwbcode==	19
					replace wbcode="BGD" if nwbcode==	20
					replace wbcode="BGR" if nwbcode==	21
					replace wbcode="BHR" if nwbcode==	22
					replace wbcode="BIH" if nwbcode==	24
					replace wbcode="BLR" if nwbcode==	25
					replace wbcode="BOL" if nwbcode==	29
					replace wbcode="BRA" if nwbcode==	30
					replace wbcode="BRB" if nwbcode==	31
					replace wbcode="BWA" if nwbcode==	36
					replace wbcode="CAF" if nwbcode==	37
					replace wbcode="CAN" if nwbcode==	38
					replace wbcode="CHE" if nwbcode==	39
					replace wbcode="CHL" if nwbcode==	40
					replace wbcode="CHN" if nwbcode==	41
					replace wbcode="CIV" if nwbcode==	42
					replace wbcode="CMR" if nwbcode==	43
					replace wbcode="COG" if nwbcode==	45
					replace wbcode="COL" if nwbcode==	46
					replace wbcode="COM" if nwbcode==	47
					replace wbcode="CPV" if nwbcode==	48
					replace wbcode="CRI" if nwbcode==	49
					replace wbcode="CUB" if nwbcode==	51
					replace wbcode="CYP" if nwbcode==	52
					replace wbcode="CZE" if nwbcode==	53
					replace wbcode="DEU" if nwbcode==	55
					replace wbcode="DJI" if nwbcode==	56
					replace wbcode="DNK" if nwbcode==	58
					replace wbcode="DOM" if nwbcode==	59
					replace wbcode="DZA" if nwbcode==	60
					replace wbcode="ECU" if nwbcode==	61
					replace wbcode="EGY" if nwbcode==	62
					replace wbcode="ESP" if nwbcode==	64
					replace wbcode="EST" if nwbcode==	65
					replace wbcode="ETH" if nwbcode==	66
					replace wbcode="FIN" if nwbcode==	67
					replace wbcode="FRA" if nwbcode==	69
					replace wbcode="GAB" if nwbcode==	71
					replace wbcode="GBR" if nwbcode==	72
					replace wbcode="GEO" if nwbcode==	73
					replace wbcode="GHA" if nwbcode==	74
					replace wbcode="GIN" if nwbcode==	76
					replace wbcode="GMB" if nwbcode==	77
					replace wbcode="GNB" if nwbcode==	78
					replace wbcode="GNQ" if nwbcode==	79
					replace wbcode="GRC" if nwbcode==	80
					replace wbcode="GTM" if nwbcode==	83
					replace wbcode="HKG" if nwbcode==	87
					replace wbcode="HND" if nwbcode==	89
					replace wbcode="HRV" if nwbcode==	91
					replace wbcode="HTI" if nwbcode==	92
					replace wbcode="HUN" if nwbcode==	93
					replace wbcode="IDN" if nwbcode==	95
					replace wbcode="IND" if nwbcode==	96
					replace wbcode="IRL" if nwbcode==	97
					replace wbcode="IRN" if nwbcode==	98
					replace wbcode="IRQ" if nwbcode==	99
					replace wbcode="ISL" if nwbcode==	100
					replace wbcode="ISR" if nwbcode==	101
					replace wbcode="ITA" if nwbcode==	102
					replace wbcode="JAM" if nwbcode==	103
					replace wbcode="JOR" if nwbcode==	104
					replace wbcode="JPN" if nwbcode==	105
					replace wbcode="KAZ" if nwbcode==	106
					replace wbcode="KEN" if nwbcode==	107
					replace wbcode="KGZ" if nwbcode==	108
					replace wbcode="KHM" if nwbcode==	109
					replace wbcode="KOR" if nwbcode==	112
					replace wbcode="KWT" if nwbcode==	113
					replace wbcode="LAO" if nwbcode==	115
					replace wbcode="LBN" if nwbcode==	116
					replace wbcode="LBR" if nwbcode==	117
					replace wbcode="LBY" if nwbcode==	118
					replace wbcode="LKA" if nwbcode==	120
					replace wbcode="LSO" if nwbcode==	121
					replace wbcode="LTU" if nwbcode==	122
					replace wbcode="LUX" if nwbcode==	123
					replace wbcode="LVA" if nwbcode==	124
					replace wbcode="MAR" if nwbcode==	126
					replace wbcode="MDA" if nwbcode==	128
					replace wbcode="MDG" if nwbcode==	129
					replace wbcode="MEX" if nwbcode==	132
					replace wbcode="MLI" if nwbcode==	135
					replace wbcode="MLT" if nwbcode==	136
					replace wbcode="MMR" if nwbcode==	137
					replace wbcode="MNE" if nwbcode==	138
					replace wbcode="MNG" if nwbcode==	139
					replace wbcode="MOZ" if nwbcode==	140
					replace wbcode="MRT" if nwbcode==	141
					replace wbcode="MUS" if nwbcode==	143
					replace wbcode="MWI" if nwbcode==	144
					replace wbcode="MYS" if nwbcode==	145
					replace wbcode="NAM" if nwbcode==	146
					replace wbcode="NER" if nwbcode==	148
					replace wbcode="NGA" if nwbcode==	149
					replace wbcode="NIC" if nwbcode==	150
					replace wbcode="NLD" if nwbcode==	151
					replace wbcode="NOR" if nwbcode==	152
					replace wbcode="NPL" if nwbcode==	153
					replace wbcode="NZL" if nwbcode==	156
					replace wbcode="OMN" if nwbcode==	158
					replace wbcode="PAK" if nwbcode==	159
					replace wbcode="PAN" if nwbcode==	160
					replace wbcode="PER" if nwbcode==	161
					replace wbcode="PHL" if nwbcode==	162
					replace wbcode="POL" if nwbcode==	165
					replace wbcode="PRK" if nwbcode==	168
					replace wbcode="PRT" if nwbcode==	170
					replace wbcode="PRY" if nwbcode==	171
					replace wbcode="QAT" if nwbcode==	176
					replace wbcode="RUS" if nwbcode==	179
					replace wbcode="RWA" if nwbcode==	180
					replace wbcode="SAU" if nwbcode==	181
					replace wbcode="SDN" if nwbcode==	183
					replace wbcode="SEN" if nwbcode==	184
					replace wbcode="SGP" if nwbcode==	185
					replace wbcode="SLE" if nwbcode==	187
					replace wbcode="SLV" if nwbcode==	188
					replace wbcode="STP" if nwbcode==	195
					replace wbcode="SVK" if nwbcode==	198
					replace wbcode="SVN" if nwbcode==	199
					replace wbcode="SWE" if nwbcode==	200
					replace wbcode="SWZ" if nwbcode==	201
					replace wbcode="SYC" if nwbcode==	203
					replace wbcode="SYR" if nwbcode==	204
					replace wbcode="TCD" if nwbcode==	205
					replace wbcode="TGO" if nwbcode==	206
					replace wbcode="THA" if nwbcode==	207
					replace wbcode="TJK" if nwbcode==	208
					replace wbcode="TKM" if nwbcode==	209
					replace wbcode="TTO" if nwbcode==	214
					replace wbcode="TUN" if nwbcode==	215
					replace wbcode="TUR" if nwbcode==	216
					replace wbcode="TZA" if nwbcode==	220
					replace wbcode="UGA" if nwbcode==	221
					replace wbcode="UKR" if nwbcode==	222
					replace wbcode="URY" if nwbcode==	223
					replace wbcode="USA" if nwbcode==	224
					replace wbcode="UZB" if nwbcode==	225
					replace wbcode="VEN" if nwbcode==	228
					replace wbcode="VNM" if nwbcode==	229
					replace wbcode="YEM" if nwbcode==	235
					replace wbcode="ZAF" if nwbcode==	238
					replace wbcode="ZMB" if nwbcode==	240
					replace wbcode="ZWE" if nwbcode==	241
					order nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
					keep nwbcode wbcode startyear endyear beta1 flips democratisations exposure total_obs nonexposure
					gen regime_def = "$democracy"
					local startyearx "$startyear"
					local endyearx "$endyear" 
					local democracyx "$democracy" 
					save "$outputfolder/PCDID_Maddison2024_Drilling_`startyearx'_`endyearx'_st_multiple_`democracyx'_diagnostics_`f'.dta", replace
					
					local i = `i'+1
					
						}
					* end of quietly
					
				di in ye _n "*****************************************************************"
				di in ye "*** Single regime cutoff defined by " in gr "dem_$democracy" 
				di in ye "*****************************************************************"
				di in ye "*** PCDID (MG) estimate with" in gr " `f' " in ye "common factor(s) included      ***"
				di in ye "*****************************************************************"
*				di in ye "Bai & NG: PCp1 " %5.3f PC1_`f' 
				di in ye "*****************************************************************"
				di in ye "Nominal sample period: " in gr "$startyear" in ye " to " in gr "$endyear "
				di in ye "resulting in " in gr "`treated'" in ye " treated countries with " ///
					in gr "`obs'" in ye " observations. "
				di in ye "MSE: " in gr %5.3f mse
				di in ye "Control sample: " in gr "`below'" in ye " countries with " in gr "`belowobs'" in ye " obs." 
				di in ye "Rob. Beta: " in gr %5.3f b_1 in ye "  SE: " in gr %5.3f se_1 ///
					in ye "  t-ratio: " in gr %5.3f b_1/se_1 in ye "  p: " in gr %5.3f 2*(normal(-abs(b_1/se_1)))
				di in ye "*****************************************************************"
				di in ye "Alpha Test: " in gr  %5.3f alphaest1 in ye " t-stat: " in gr  %5.3f alphat1
				di in ye "*****************************************************************"
				di in ye "*** Test of controls " in gr %5.2f `chi2_stat' in ye " (" in gr %5.3f `chi2_p' in ye ").      ***"
				di in ye "*****************************************************************"
				local ii = `i'-1
				di in red "`ii'"
				
				restore	
				
			}
			* end of the factor loop
		}
		* end of dem indicator loop
	}
	* end of start year loop
}
* end of end year loop

local mytitle "PCDID"
local esttabformat  "nodepvars lines compress label nogaps numbers star(* 0.10 ** 0.05 *** 0.01)"
local esttabstats "stats(start_year final_year fac treated obs controlN controlobs x_test_chi2 x_test_p alpha_b1 alpha_t1, fmt(0 0 0 0 0 0 0 2 2 3 3))"
local esttabkeep "keep($demo_list2 $xvars2) order($demo_list2 $xvars2)"
 esttab model* using "$outputfolder/20240605_PCDID_1959_2018_Maddison_Drilling_High_Mid_level_multiple_Diagnostics.csv", b(3) se(3) abs br ///
 	`esttabformat' `esttabkeep' `esttabstats' style(tab) replace nodep

exit	
